In [ ]:

    

%matplotlib inline
from math import *
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from QuantLib import *
from PyFin.Math.Distributions import CumulativeNormalDistribution

plt.style.use('fivethirtyeight')


    

In [ ]:

    

def bs_theoretical_price(payoff, spot, ttm, volatility, rf_rate, finance_rate=None):
    
    if not finance_rate:
        finance_rate = rf_rate
    
    forward = spot * exp(finance_rate * ttm)
    std_dev = volatility * sqrt(ttm)
    discount = exp(-rf_rate * ttm)
    return blackFormula(payoff.optionType(), payoff.strike(), forward, std_dev, discount)


    

In [ ]:

    

_dist = CumulativeNormalDistribution()

def _exercise(payoff, spot):
    return payoff(spot)

_exercise = np.frompyfunc(_exercise, 2, 1)


def _create_path(rsg, ln_spot, drift, diffusion, delta_t):
    rnd = np.array(rsg.nextSequence().value())
    in_c = delta_t * drift + np.sqrt(delta_t) * diffusion * rnd
    inc_c_cum = np.cumsum(in_c)
    return np.exp(np.concatenate(([ln_spot], ln_spot + inc_c_cum)))


def _bs_delta(option_type, finance_rate, volatility, ttm, spot, strike):
    money_ness = log(spot / strike)
    drift = (finance_rate + 0.5 * (volatility ** 2)) * ttm
    d1 = (money_ness + drift) / volatility / sqrt(ttm)
    call_delta =  _dist(d1)
    
    if option_type == Option.Call:
        return call_delta
    elif option_type == Option.Put:
        return call_delta - 1.

_bs_delta = np.frompyfunc(_bs_delta, 6, 1)


def _hedging_on_path(payoff, ttm, time_grids, spot_path, volatility, inflations, rf_rate, finance_rate, trading_cost):
    delta_t = time_grids[0] - time_grids[1]
    deltas = _bs_delta(payoff.optionType(), finance_rate, volatility, time_grids, spot_path[:-1], payoff.strike())
    borrows = spot_path[:-1] * deltas
    finance_cost = borrows * finance_rate * delta_t
    stock_pnl = deltas * (spot_path[1:] - spot_path[:-1])
    trading_slipge = np.abs(np.concatenate(([deltas[0]], np.diff(deltas)))) * spot_path[:-1] * trading_cost
    hedging_pnl = ((stock_pnl - finance_cost - trading_slipge) * inflations).sum()
    
    exercise_pnl = payoff(spot_path[-1])
    total_cost = hedging_pnl - exercise_pnl
    
    return exp(-rf_rate * ttm) * total_cost


class HedgeAnalysor(object):
    
    def __init__(self,
                 payoff,
                 trading_cost=0.):
        self.payoff = payoff
        self.trading_cost = trading_cost
    
    @staticmethod
    def _prepare_parameters(rf_rate, finance_rate, underlying_risk_return):
        if finance_rate is None:
            finance_rate = rf_rate
        
        if underlying_risk_return is None:
            underlying_risk_return = rf_rate
        
        return finance_rate, underlying_risk_return

    def exercise(self, spots):
        return _exercise(self.payoff, spots)
    
    def _hedge_path(self, rsg, ttm, time_grids, ln_spot, volatility, drift, diffusion, delta_t, inflations, rf_rate, finance_rate):
        spot_path = _create_path(rsg, ln_spot, drift, diffusion, delta_t)
        return _hedging_on_path(self.payoff, ttm, time_grids, spot_path, volatility, inflations, rf_rate, finance_rate, self.trading_cost)
        
    
    def hedge_cost(self,
                   rf_rate,
                   volatility,
                   realized_vol,
                   ttm,
                   finance_rate=None,
                   underlying_risk_return=None,
                   spot=1.,
                   time_steps=50,
                   simulations=100,
                   seed=20):
        rng = MersenneTwisterUniformRng(seed)
        rsg = MersenneTwisterUniformRsg(dimensionality=time_steps,
                                        rng=rng)
        rsg = InvCumulativeMersenneTwisterGaussianRsg(rsg)
        
        finance_rate, underlying_risk_return = self._prepare_parameters(rf_rate,
                                                                        finance_rate,
                                                                        underlying_risk_return)
        print("risk free: {0:.02f}%".format(rf_rate*100))
        print('finance rate: {0:.02f}%'.format(finance_rate*100))
        print('underlying risk return: {0:.02f}%'.format(underlying_risk_return*100))
        print('implied vol: {0:.02f}%'.format(volatility))
        print('realized vol: {0:.02f}%'.format(realized_vol))
        print('number of simulations: {0}'.format(simulations))
        print('time to maturity: {0} yrs.'.format(ttm))
        print('time steps: {0}'.format(time_steps))
        print('trading cost: {0:0.4f}%'.format(self.trading_cost*100))
        print('payoff: type={0}, k={1}'.format(self.payoff.optionType(), self.payoff.strike()))
        
        ln_spot = log(spot)
        delta_t = ttm / time_steps

        drift = underlying_risk_return - 0.5 * (realized_vol ** 2)
        diffusion = realized_vol
        time_grids = np.linspace(ttm, 0, num=time_steps, endpoint=False)
        inflations = np.exp(finance_rate * (time_grids - delta_t))
        
        hedging_cost_batch = np.zeros(simulations)
        
        for i in range(simulations):
            hedging_cost = self._hedge_path(rsg,
                                            ttm,
                                            time_grids,
                                            ln_spot,
                                            volatility,
                                            drift,
                                            diffusion,
                                            delta_t,
                                            inflations,
                                            rf_rate,
                                            finance_rate)
            hedging_cost_batch[i] = hedging_cost
        
        return -hedging_cost_batch


    

In [ ]:

    

rf_rate = 0.04
finance_rate = 0.06
underlying_risk_return = -0.15
strike = 1.
volatility = 0.30
realized_vol = 0.30
ttm = 0.25
spot = 1.
trading_cost = 0.0
simulations = 50000
time_steps = int(ttm * 250)

payoff = PlainVanillaPayoff(Option.Call, strike)
hf = HedgeAnalysor(payoff, trading_cost)


    

In [ ]:

    

%%time

res = hf.hedge_cost(rf_rate=rf_rate,
                    volatility=volatility,
                    realized_vol=realized_vol,
                    ttm=ttm,
                    spot=spot,
                    finance_rate=finance_rate,
                    underlying_risk_return=underlying_risk_return,
                    simulations=simulations,
                    time_steps=time_steps)
bs_price = bs_theoretical_price(payoff, spot, ttm, volatility, rf_rate, finance_rate)


    

In [ ]:

    

fig, axes = plt.subplots(1, 1, figsize=(12, 6), sharex=True)
axes = [axes]

for i, ax in enumerate(axes):
    ax.hist(res, bins=50)
    ax.axvline(x=bs_price, color='red', linestyle='dashed', label='Theoretical: {0:.02f}%'.format(bs_price*100))
    ax.axvline(x=res.mean(), color='green', linestyle='dashed', label='Hedging Mean')
    ax.axvline(x=np.percentile(res, 1), color='yellow', linestyle='dashed', label='Hedging per. 1%')
    ax.axvline(x=np.percentile(res, 5), color='black', linestyle='dashed', label='Hedging per. 5%')
    ax.axvline(x=np.percentile(res, 95), color='black', linestyle='dashed', label='Hedging per. 95%')
    ax.axvline(x=np.percentile(res, 99), color='yellow', linestyle='dashed', label='Hedging per. 99%')
    ax.set_title("Hedging v.s. Theoretical (tc = {0}%, r_vol = {1:0.2f}%, k={2})".format(trading_cost*100,
                                                                                                    realized_vol*100,
                                                                                                strike))
    ax.legend()


    

In [ ]:

    

# produce the table
index = ['理论', '对冲（平均）', '对冲（分位数 1%）', '对冲（分位数 5%）', '对冲（分位数 95%）', '对冲（分位数 99%）']
values = np.array([bs_price, res.mean(), np.percentile(res, 1), np.percentile(res, 5), np.percentile(res, 95), np.percentile(res, 99)])
rel_values = values / values[0] - 1.
df = pd.DataFrame(data={'成本': values, '相对': rel_values}, index=index)
df


    

In [ ]:

    

rf_rate = rf_rate
finance_rate = finance_rate
underlying_risk_return = underlying_risk_return
strike = 1.
volatility = 0.30
realized_vol = 0.30
ttm = 0.25
spot = 1.
trading_cost = 0.0015
simulations = 50000
time_steps = int(ttm * 250)

payoff = PlainVanillaPayoff(Option.Call, strike)
hf = HedgeAnalysor(payoff, trading_cost=trading_cost)


    

In [ ]:

    

strike_scenarios =[0.9, 1.0, 1.1]


    

In [ ]:

    

simulations_res = []
bs_prices_res = []

for i, this_strike in enumerate(strike_scenarios):
    print("\nScenarios {0} ......".format(i+1))
    this_payoff = PlainVanillaPayoff(Option.Call, this_strike)
    this_bs_price = bs_theoretical_price(this_payoff, spot, ttm, volatility, rf_rate, finance_rate)
    hf = HedgeAnalysor(this_payoff, trading_cost)
    path_res = hf.hedge_cost(rf_rate=rf_rate,
                             volatility=volatility,
                             realized_vol=realized_vol,
                             time_steps=time_steps,
                             ttm=ttm,
                             finance_rate=finance_rate,
                             underlying_risk_return=underlying_risk_return,
                             spot=spot,
                             simulations=simulations)
    
    simulations_res.append(path_res)
    bs_prices_res.append(this_bs_price)


    

In [ ]:

    

fig, axes = plt.subplots(len(strike_scenarios), 1, figsize=(12, 6 * len(strike_scenarios)), sharex=True)

if not hasattr(axes, '__iter__'):
    axes = [axes]

for i, ax in enumerate(axes):
    res = simulations_res[i]
    ax.hist(res, bins=50)
    ax.axvline(x=bs_prices_res[i], color='red', linestyle='dashed', label='Theoretical: {0:.02f}%'.format(bs_prices_res[i]))
    ax.axvline(x=res.mean(), color='green', linestyle='dashed', label='Hedging Mean')
    ax.axvline(x=np.percentile(res, 1), color='yellow', linestyle='dashed', label='Hedging per. 1%')
    ax.axvline(x=np.percentile(res, 5), color='black', linestyle='dashed', label='Hedging per. 5%')
    ax.axvline(x=np.percentile(res, 95), color='black', linestyle='dashed', label='Hedging per. 95%')
    ax.axvline(x=np.percentile(res, 99), color='yellow', linestyle='dashed', label='Hedging per. 99%')
    ax.set_title("Hedging v.s. Theoretical (tc = {0:.2f}%, r_vol = {1:0.2f}%, k={2})".format(trading_cost*100,
                                                                                                    realized_vol*100,
                                                                                                       strike_scenarios[i]))
    ax.legend()


    

In [ ]:

    

# produce the table
index = ['理论', '对冲（平均）', '对冲（分位数 1%）', '对冲（分位数 5%）', '对冲（分位数 95%）', '对冲（分位数 99%）']
col_names = ['{0:.0f}%'.format(t*100) for t in strike_scenarios]

values = np.zeros((len(index), len(col_names)))
for j, res in enumerate(simulations_res):
    this_values = np.array([bs_prices_res[j], res.mean(), np.percentile(res, 1), np.percentile(res, 5), np.percentile(res, 95), np.percentile(res, 99)])
    values[:, j] = this_values
    
df = pd.DataFrame(data=values, columns=col_names, index=index)
df


    

In [ ]:

    

trading_cost_scenarios =[0.0010, 0.0015, 0.0020]


    

In [ ]:

    

simulations_res = []
bs_price = bs_theoretical_price(payoff, spot, ttm, volatility, rf_rate, finance_rate)

for i, this_trading_cost in enumerate(trading_cost_scenarios):
    print("\nScenarios {0} ......".format(i+1))
    hf = HedgeAnalysor(payoff, this_trading_cost)
    path_res = hf.hedge_cost(rf_rate=rf_rate,
                             volatility=volatility,
                             realized_vol=realized_vol,
                             time_steps=time_steps,
                             ttm=ttm,
                             finance_rate=finance_rate,
                             underlying_risk_return=underlying_risk_return,
                             spot=spot,
                             simulations=simulations)
    
    simulations_res.append(path_res)


    

In [ ]:

    

fig, axes = plt.subplots(len(trading_cost_scenarios), 1, figsize=(12, 6 * len(trading_cost_scenarios)), sharex=True)

if not hasattr(axes, '__iter__'):
    axes = [axes]

for i, ax in enumerate(axes):
    res = simulations_res[i]
    ax.hist(res, bins=50)
    ax.axvline(x=bs_price, color='red', linestyle='dashed', label='Theoretical: {0:.02f}%'.format(bs_price*100))
    ax.axvline(x=res.mean(), color='green', linestyle='dashed', label='Hedging Mean')
    ax.axvline(x=np.percentile(res, 1), color='yellow', linestyle='dashed', label='Hedging per. 1%')
    ax.axvline(x=np.percentile(res, 5), color='black', linestyle='dashed', label='Hedging per. 5%')
    ax.axvline(x=np.percentile(res, 95), color='black', linestyle='dashed', label='Hedging per. 95%')
    ax.axvline(x=np.percentile(res, 99), color='yellow', linestyle='dashed', label='Hedging per. 99%')
    ax.set_title("Hedging v.s. Theoretical (tc = {0:.2f}%, r_vol = {1:0.2f}%, k={2})".format(trading_cost_scenarios[i]*100,
                                                                                                    realized_vol*100,
                                                                                                       strike))
    ax.legend()


    

In [ ]:

    

# produce the table
index = ['理论', '对冲（平均）', '对冲（分位数 1%）', '对冲（分位数 5%）', '对冲（分位数 95%）', '对冲（分位数 99%）']
col_names = ['{0:.2f}%'.format(t*100) for t in trading_cost_scenarios]

values = np.zeros((len(index), len(col_names)))
for j, res in enumerate(simulations_res):
    this_values = np.array([bs_price, res.mean(), np.percentile(res, 1), np.percentile(res, 5), np.percentile(res, 95), np.percentile(res, 99)])
    values[:, j] = this_values
    
df = pd.DataFrame(data=values, columns=col_names, index=index)
df


    

In [ ]:

    

volatility_scenarios = [0.20, 0.30, 0.40]


    

In [ ]:

    

simulations_res = []
bs_price = bs_theoretical_price(payoff, spot, ttm, volatility, rf_rate, finance_rate)
hf = HedgeAnalysor(payoff, trading_cost)

for i, this_volatility in enumerate(volatility_scenarios):
    print("\nScenarios {0} ......".format(i+1))
    path_res = hf.hedge_cost(rf_rate=rf_rate,
                             volatility=volatility,
                             realized_vol=this_volatility,
                             time_steps=time_steps,
                             ttm=ttm,
                             finance_rate=finance_rate,
                             underlying_risk_return=underlying_risk_return,
                             spot=spot,
                             simulations=simulations)
    
    simulations_res.append(path_res)


    

In [ ]:

    

fig, axes = plt.subplots(len(volatility_scenarios), 1, figsize=(12, 6 * len(volatility_scenarios)), sharex=True)

if not hasattr(axes, '__iter__'):
    axes = [axes]

for i, ax in enumerate(axes):
    res = simulations_res[i]
    ax.hist(res, bins=50)
    ax.axvline(x=bs_price, color='red', linestyle='dashed', label='Theoretica: {0:.02f}%'.format(bs_price*100))
    ax.axvline(x=res.mean(), color='green', linestyle='dashed', label='Hedging Mean')
    ax.axvline(x=np.percentile(res, 1), color='yellow', linestyle='dashed', label='Hedging per. 1%')
    ax.axvline(x=np.percentile(res, 5), color='black', linestyle='dashed', label='Hedging per. 5%')
    ax.axvline(x=np.percentile(res, 95), color='black', linestyle='dashed', label='Hedging per. 95%')
    ax.axvline(x=np.percentile(res, 99), color='yellow', linestyle='dashed', label='Hedging per. 99%')
    ax.set_title("Hedging v.s. Theoretical (tc = {0:.2f}%, r_vol = {1:0.2f}%, k={2})".format(trading_cost*100,
                                                                                                    volatility_scenarios[i]*100,
                                                                                                       strike))
    ax.legend()


    

In [ ]:

    

# produce the table
index = ['理论', '对冲（平均）', '对冲（分位数 1%）', '对冲（分位数 5%）', '对冲（分位数 95%）', '对冲（分位数 99%）']
col_names = ['{0:.2f}%'.format(v*100) for v in volatility_scenarios]

values = np.zeros((len(index), len(col_names)))
for j, res in enumerate(simulations_res):
    this_values = np.array([bs_price, res.mean(), np.percentile(res, 1), np.percentile(res, 5), np.percentile(res, 95), np.percentile(res, 99)])
    values[:, j] = this_values
    
df = pd.DataFrame(data=values, columns=col_names, index=index)
df


    

In [ ]:

    




    

In [ ]: